## Government polling in times of crises: When capacity meets incentives
## JEPP
## Tinette Schnatterer & Anja Durovic
## 20250701

library(mlogit)
library(texreg)

ReplicationCrises_GermanGovernmentPolls <- read.csv("ReplicationCrises_GermanGovernmentPolls.csv", header= T, sep= ",")

#### Main Model: Conditional Logit Regression Without FE------------- 

ReplicationCrises_GermanGovernmentPolls <- mlogit.data(ReplicationCrises_GermanGovernmentPolls, 
                                                       choice = "Choice",   # Binary outcome: 1 if issue was chosen
                                                       shape = "long", 
                                                       alt.var = "PolicyIssue",  # Alternative-specific variable
                                                       chid.var = "id") # Choice-set identifier


ReplicationCrises_GermanGovernmentPolls$LegResponsibility <- as.factor(ReplicationCrises_GermanGovernmentPolls$LegResponsibility)

model1CL <- mlogit(Choice ~  LegResponsibility +  MediaSalience + GovernmentPriorities + IssueOwnership | 0, 
                   data = ReplicationCrises_GermanGovernmentPolls)
summary(model1CL)

model2CL <- mlogit(Choice ~  LegResponsibility + PolicyType +  MediaSalience + GovernmentPriorities + IssueOwnership| 0, 
                   data = ReplicationCrises_GermanGovernmentPolls)
summary(model2CL)

model3CL <- mlogit(Choice ~  LegResponsibility + PolicyType + IssueSalience +  MediaSalience + GovernmentPriorities + IssueOwnership | 0, 
                   data = ReplicationCrises_GermanGovernmentPolls)
summary(model3CL)

model4CL <-  mlogit(Choice  ~ IssueSalience  * LegResponsibility +LegResponsibility + PolicyType + IssueSalience + GovernmentPriorities + MediaSalience + IssueOwnership| 0, 
                    data = ReplicationCrises_GermanGovernmentPolls)
summary(model4CL)

model5CL <-  mlogit(Choice  ~ IssueSalience  * PolicyType + LegResponsibility + PolicyType + IssueSalience + GovernmentPriorities+ MediaSalience + IssueOwnership | 0, 
                    data = ReplicationCrises_GermanGovernmentPolls)
summary(model5CL)
exp(coef(model5CL))

texreg(list(model1CL, model2CL, model3CL, model4CL, model5CL), 
       file = "modelsCL_summary.tex", 
       caption = "Issue selection in German government polls (conditional logit regression", 
       label = "tab:models", 
       use.packages = FALSE)



get_transformed_output <- function(model) {
  coefs <- coef(model)
  ses <- sqrt(diag(vcov(model)))
  
  # Exponentiated coefficients (Odds Ratios)
  ORs <- exp(coefs)
  
  # Delta method: transformed standard errors
  OR_ses <- ORs * ses
  
  # Z-values and corresponding p-values based on original coefficients (log-odds)
  zvals <- coefs / ses
  pvals <- 2 * (1 - pnorm(abs(zvals)))
  
  return(list(coef = ORs, se = OR_ses, pval = pvals))
}


# Apply to each model
models <- list(model1CL, model2CL, model3CL, model4CL, model5CL)
transformed <- lapply(models, get_transformed_output)

texreg(models,
       file = "multiple_modelsCLOdds.tex",
       label = "tab:multiple_modelsCL",
       use.packages = FALSE,
       dcolumn = TRUE,
       stars = c(0.1, 0.05, 0.01),
       custom.model.names = c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5"),
       override.coef = lapply(transformed, function(x) x$coef),
       override.se   = lapply(transformed, function(x) x$se),
       override.pval = lapply(transformed, function(x) x$pval)
)




texreg(list(model1CL, model2CL, model3CL, model4CL, model5CL), 
       file = "multiple_modelsOddsCL2.tex",
       #caption = "Comparison of Multiple ModelsOdds",
       label = "tab:multiple_models",
       use.packages = FALSE,
       dcolumn = TRUE, 
       stars = c(0.1, 0.05, 0.01),
       custom.model.names = c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5"),
       transform = exp,  # <-- this gives you odds ratios!
       override.se = lapply(list(model1CL, model2CL, model3CL, model4CL, model5CL), function(m) sqrt(diag(vcov(m))))
)



#### Mixed logit models----

### Lagrange Test-----
names(ReplicationCrises_GermanGovernmentPolls)

mixed_logit_modelT <- mlogit(Choice ~  LegResponsibility + PolicyType + IssueSalience + MediaSalience + GovernmentPriorities  +IssueOwnership | 0, 
                             data = ReplicationCrises_GermanGovernmentPolls)

summary(mixed_logit_modelT)
------
  
  mixed_logit_modelTR1 <- mlogit(Choice ~LegResponsibility + PolicyType + IssueSalience + MediaSalience + GovernmentPriorities  +IssueOwnership| 0, 
                                 data = ReplicationCrises_GermanGovernmentPolls,rpar = c(IssueSalience= "n"), R = 100)

summary(mixed_logit_modelTR1)

mixed_logit_modelTR2 <- mlogit(Choice ~  LegResponsibility + PolicyType + MediaSalience + GovernmentPriorities  +IssueOwnership | 0, 
                               data = ReplicationCrises_GermanGovernmentPolls,rpar = c(MediaSalience= "n"), R = 100)

summary(mixed_logit_modelTR2)

mixed_logit_modelTR3 <- mlogit(Choice ~  LegResponsibility + PolicyType + IssueSalience + MediaSalience + GovernmentPriorities  +IssueOwnership | 0, 
                               data = ReplicationCrises_GermanGovernmentPolls,rpar = c(GovernmentPriorities= "n"), R = 100)

summary(mixed_logit_modelTR3)


mixed_logit_modelTR4 <- mlogit(Choice ~  LegResponsibility + PolicyType +IssueSalience + MediaSalience + GovernmentPriorities  + IssueOwnership | 0, 
                               data = ReplicationCrises_GermanGovernmentPolls,rpar = c(GovernmentPriorities= "n", IssueSalience = "n", 
                                                                                       MediaSalience = "n"), R = 100)

summary(mixed_logit_modelTR4)

mixed_logit_modelTR5 <- mlogit(Choice ~  LegResponsibility*IssueSalience + LegResponsibility + PolicyType + IssueSalience + MediaSalience + GovernmentPriorities + IssueOwnership| 0, 
                               data = ReplicationCrises_GermanGovernmentPolls,rpar = c(GovernmentPriorities= "n", IssueSalience = "n", 
                                                                                       MediaSalience = "n"), R = 100)

summary(mixed_logit_modelTR5)

mixed_logit_modelTR6 <- mlogit(Choice ~  PolicyType*IssueSalience + LegResponsibility + PolicyType + IssueSalience + MediaSalience + GovernmentPriorities + IssueOwnership | 0, 
                               data = ReplicationCrises_GermanGovernmentPolls,rpar = c(GovernmentPriorities= "n", IssueSalience = "n", 
                                                                                       MediaSalience = "n"), R = 100)

summary(mixed_logit_modelTR6)

lrt_stat <- -2 * (logLik(mixed_logit_modelT ) - logLik(mixed_logit_modelTR1))
lrt_df <- length(coef(mixed_logit_modelTR1)) - length(coef(mixed_logit_modelT ))
lrt_p <- pchisq(lrt_stat, df = lrt_df, lower.tail = FALSE)

# Print the results
cat("Likelihood Ratio Test Statistic:", lrt_stat, "\n")
cat("Degrees of Freedom:", lrt_df, "\n")
cat("p-value:", lrt_p, "\n")

lrt_stat <- -2 * (logLik(mixed_logit_modelT ) - logLik(mixed_logit_modelTR2))
lrt_df <- length(coef(mixed_logit_modelTR2)) - length(coef(mixed_logit_modelT ))
lrt_p <- pchisq(lrt_stat, df = lrt_df, lower.tail = FALSE)

# Print the results
cat("Likelihood Ratio Test Statistic:", lrt_stat, "\n")
cat("Degrees of Freedom:", lrt_df, "\n")
cat("p-value:", lrt_p, "\n")



lrt_stat <- -2 * (logLik(mixed_logit_modelT ) - logLik(mixed_logit_modelTR3))
lrt_df <- length(coef(mixed_logit_modelTR3)) - length(coef(mixed_logit_modelT ))
lrt_p <- pchisq(lrt_stat, df = lrt_df, lower.tail = FALSE)

# Print the results
cat("Likelihood Ratio Test Statistic:", lrt_stat, "\n")
cat("Degrees of Freedom:", lrt_df, "\n")
cat("p-value:", lrt_p, "\n")

lrt_stat <- -2 * (logLik(mixed_logit_modelT ) - logLik(mixed_logit_modelTR4))
lrt_df <- length(coef(mixed_logit_modelTR4)) - length(coef(mixed_logit_modelT ))
lrt_p <- pchisq(lrt_stat, df = lrt_df, lower.tail = FALSE)

# Print the results
cat("Likelihood Ratio Test Statistic:", lrt_stat, "\n")
cat("Degrees of Freedom:", lrt_df, "\n")
cat("p-value:", lrt_p, "\n")

lrt_stat <- -2 * (logLik(mixed_logit_modelT ) - logLik(mixed_logit_modelTR5))
lrt_df <- length(coef(mixed_logit_modelTR5)) - length(coef(mixed_logit_modelT ))
lrt_p <- pchisq(lrt_stat, df = lrt_df, lower.tail = FALSE)

# Print the results
cat("Likelihood Ratio Test Statistic:", lrt_stat, "\n")
cat("Degrees of Freedom:", lrt_df, "\n")
cat("p-value:", lrt_p, "\n")

lrt_stat <- -2 * (logLik(mixed_logit_modelT ) - logLik(mixed_logit_modelTR6))
lrt_df <- length(coef(mixed_logit_modelTR6)) - length(coef(mixed_logit_modelT ))
lrt_p <- pchisq(lrt_stat, df = lrt_df, lower.tail = FALSE)

# Print the results
cat("Likelihood Ratio Test Statistic:", lrt_stat, "\n")
cat("Degrees of Freedom:", lrt_df, "\n")
cat("p-value:", lrt_p, "\n")


texreg(list(mixed_logit_modelTR4, mixed_logit_modelTR5, mixed_logit_modelTR6), 
       file = "mixed_logit.tex", 
       caption = "Mixed Logit Models with Random Effects", 
       label = "tab:models", 
       use.packages = FALSE)
